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Abstract. Primordial non-Gaussianity can lead to a scale-dependent bias in the density of collapsed 
halos relative to the underlying matter density. The galaxy power spectrum already provides con¬ 
straints on local-type primordial non-Gaussianity complementary those from the cosmic microwave 
background (CMB), while the bispectrum contains additional shape information and has the potential 
to outperform CMB constraints in future. We develop the bias model for the halo density contrast 
in the presence of local-type primordial non-Gaussianity, deriving a bivariate expansion up to sec¬ 
ond order in terms of the local linear matter density contrast and the local gravitational potential 
in Lagrangian coordinates. Nonlinear evolution of the matter density introduces a non-local tidal 
term in the halo model. Furthermore, the presence of local-type non-Gaussianity in the Lagrangian 
frame leads to a novel non-local convective term in the Eulerian frame, that is proportional to the 
displacement field when going beyond the spherical collapse approximation. We use an extended 
Press-Schechter approach to evaluate the halo mass function and thus the halo bispectrum. We show 
that including these non-local terms in the halo bispectra can lead to corrections of up to 25% for 
some configurations, on large scales or at high redshift. 
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1 Introduction 

1.1 General motivations 

Inflation is a powerful mechanism to solve the puzzle of initial conditions for the Hot Big Bang cos¬ 
mology and explain the origin of structure in our Universe. However, the exact mechanism that drove 
inflation is still a matter of debate. Different models make different predictions for the statistics of 
the primordial curvature perturbations (q n f, which can be quantified by computing the n-point cor¬ 
relation functions (Cinf(ki)... Cmf(k n ))- All models predict some departure from Gaussianity, giving 
non-vanishing correlation functions for n > 2. Since the amplitude, shape- and scale-dependence are 
model dependent, the intrinsic non-Gaussianity of primordial curvature perturbations is a powerful 
tool to discriminate among different models [1], For example, local-type non Gaussianity is a dis¬ 
tinctive signature of multi-field models [2], while equilateral and orthogonal non-Gaussianity probe 
non-slow roll dynamics through the self-interactions of the fields [3]. Thus observational features of 
the primordial curvature perturbation offer a window onto the nature of inflation [4]. 

At the time of this publication, the most accurate constraints on primordial non-Gaussianity 
(PNG) come from the cosmic microwave background (CMB). WMAP9 set the constraint on local /nl 
to be —3 < /nl < 77 at 95% CL [5], while the most recent data from the Planck satellite experiment 
requires —9.2 < / NL < 10.8 [6]. This is compatible with the predictions of simple slowly-rolling 
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single-field models which predict |/nl| < 1- The CMB is a 2D map at redshift 2 ss 1100 whose 
temperature fluctuations have been measured with high accuracy. The error bars could be further 
decreased by including the higher-resolution polarization data [7], but we are close to having fully 
exploited the constraining power of the CMB, being limited on small scales by the Silk damping and 
on large scales by the cosmic variance. The question is then which other observables may be suitable 
to obtain bounds on A/nl ~ 1. 

The large-scale structure (LSS) of the Universe provides a 3D map of the Universe, with accessible 
modes potentially going from the horizon scale ~ 10 4 Mpc down to the non-linear scale ~ 10 Mpc. 
The pioneering paper by Dalai et al. [8] (see also [9]) showed that the LSS power spectrum also 
offers a strong constraint on primordial non-Gaussianity. The idea comes from the mode coupling 
in local-type models; the short wavelength modes which drive the collapse of matter into halos are 
modulated by the long wavelength modes, effectively changing the short-mode variance from patch 
to patch of the sky. This introduces a scale-dependent correction in the bias model, which describes 
how the number of halos relates to the matter density. Since the correction is proportional to /nl/& 2 , 
the clustering of structures on large scales has the potential to discriminate among different models 
of inflation. 

It is important to emphasise that LSS constraints are independent of those from the CMB and 
apply on different scales, hence constraining also scale-dependent non-Gaussianity. This new perspec¬ 
tive on the subject has been studied in depth in a number of works [9-15]. However, constraining /nl 
from real data requires detailed understanding of the systematic errors; indeed, some of them mimic 
the PNG excess of power on large scales (see for instance [16-18]). Different techniques have been 
proposed to handle these errors [17, 19-21], as well as methods to reduce the effect of cosmic variance 
[22-27] and the statistical uncertainties [28, 29]. Constraints comparable with the WMAP bounds 
have already been achieved [10, 19-21, 30-34] and even more stringent results will come thanks to 
the next generation of experiments, like Euclid and SKA [27, 35-37], despite the fact that no survey 
has been optimized for PGN so far [38]. The novel technique introduced in [39] may optimistically 
provide A/ NL ~ 1. 

However the 2-point function is not the natural statistic in which to look for PNG, as the full 
shape information is available only in higher-order statistics, and moreover local-type models including 
higher-order non-Gaussian parameters like 3 nl are approximately degenerate with /nl in the power 
spectrum [40]. These issues naturally drive us towards studying higher-order n-points correlation 
functions, in particular the 3-point function and its Fourier transform, the bispectrum 1 . Many more 
3D triangle configurations compared to 2D CMB ones are available in fc-space, suggesting potentially 
a stronger constraining power with respect to the power spectrum. Indeed, A/nl ~ 1 is expected 
from the bispectrum [42-44], although this forecast is based on idealized surveys and ignores redshift 
space distortions (RSD). In addition, the degeneracy between /nl and $nl is broken by combining 
both the power spectrum and bispectrum constraints [45, 46]. 

There are a number of complications, which explains why few measurements of the 3-point 
function or bispectrum from real data have been obtained so far [42, 47-55]. For example, redshift- 
space distortions and the complicated mask geometry intrinsic to any survey are hard to model. 
Further, non-linearities in halo populations and the non-linear nature of Einsteins field equations 
produce non-Gaussian features in the evolved matter field that may be hard to disentangle from the 
PNG signal. 

1.2 What we do in this work 

The ability to go beyond current constraints on PNG, crucially depends on our understanding of 
all possible effects that can influence LSS measurements and how these depend on the primordial 
fluctuations. In this paper we make improvements to the bias model which relates the halo density to 
the underlying matter density, particularly focussing on an accurate description of the second-order, 
non-local and non-Gaussian effects. 

A common procedure to describe how dark matter halos trace the matter density is to assume a 
local relation between the final number density of objects with mass M at redshift z and the evolved 

1 In [41] the possibility to use the higher-order moments of LSS to constraint PNG has been explored. Although 
these avoid the complexities of the full n-point statistics, they may not be able to detect small PNG. 
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density field, known as local Eulerian biasing [56], 

oo 7 E / \ 

#h( x , A M ) = z)Y , (1.1) 

7=1 -7 ' 

where bj are the bias coefficients and is the fully non-linear density contrast in Eulerian space, 
smoothed on the mass scale M. In this picture, the halos are drawn on top of the density field 
without any memory of their past history, with possible application also to the case of PNG [42, 43, 45]. 
However, many recent results show that this biasing model is not sufficiently accurate when compared 
to simulations [57-61]. In particular, other physically motivated contributions consistent with the 
symmetries of the dark matter equations of motion [62], like a tidal term, should be included. These 
new terms break the assumption of locality implicit in eq. (1.1). Either a local or non-local Eulerian 
model offer only a parametrization of the bias, and the bias coefficients need to be fitted against data 
or N-body simulations. 

A promising alternative is to assume a physically reasonable local relation between the initial 
number density of objects and the initial density field [63] , 

00 b L (z) 

<$(q. a M ) = ~r~ (^n,M(q> z )Y > ( L2 ) 

7=1 

where S^ n M is the linearly evolving density contrast smoothed on scale M. Note that the expansion 
is performed in Lagrangian space, with the initial spatial coordinate q being related to the evolved 
Eulerian coordinate, x, through the relation 

x(q,r) = q+^(q,r) , (1.3) 

where M/, the displacement field, is the key dynamical quantity in the Lagrangian picture. Equa¬ 
tion (1.2) is known as a local Lagrangian biasing scheme and in effect assumes that the formation 
sites of halos can be identified from the initial density field 2 . The dynamics of halos is then captured 
in the transformation to the Eulerian space, by applying eq. (1.3). This local Lagrangian biasing 
scheme comes with a prescription to calculate the bias, if we know the halo mass function, and it also 
automatically leads to the presence of a tidal term when written in terms of the non-linearly evolved 
density field [57, 66]. We will follow this approach, carefully evaluating the effects of the displacement 
field on the physical quantities we will investigate. 

In presence of PNG of local type, the usual local Lagrangian bias model needs to be extended 
to account for the correlations in the initial density field. In the following we start re-deriving the 
bivariate model of [67], which remarkably shows good fits against simulations even when applied to 
the bispectrum [44, 68]. We then show that a novel non-local term arises, contributing to the number 
density of halos, which was previously neglected under the assumption of spherical collapse. We will 
quantitatively investigate how much the tidal term, and our new non-local contribution, affect the 
tree-level halo bispectrum in comparison with the reference model of [44] (hereafter BSS) 3 . 

This paper is organized as follows. In section 2, we briefly introduce the main concepts and 
notation. In particular, in section 2.1 we review the evolution of the matter density field up to second 
order and indentify a new convective term in presence of PNG, while in section 2.2 we describe the 
peak-background split. This constitutes the basis of our local Lagrangian biasing scheme, which 
will be presented in section 3. As we work in a Lagrangian approach, our result then needs to be 
transformed into the Eulerian frame in section 4. In section 5 expressions for the matter and halo 
bispectra are given, and a numerical analysis of the halo bispectrum is presented in section 6, where 
we quantify the effect of new non-local terms in LSS bispectra. We present our final conclusions 
together with a discussion of potential future developments of our work in section 7. 

2 Generalizations of this picture exists as well, see for instance [64, 65]. 

3 Alternative frameworks to include PNG of general type in a Lagrangian description exist as well. See for instance 
[69] and [70], which is based on the integrated perturbation theory developed in [65] 
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2 Basic framework 


2.1 Matter density field 

In this subsection we define the quantities and introduce the notation that we will use in what follows. 
The linearly growing mode of the density field has a simple separable form in Lagrangian space 


<Wq> z) = c(q )D(z ), 


( 2 . 1 ) 


where D{z) is the linear growth factor with D(z ) oc (l + z)^ in the matter-dominated era, normalized 
such that D( 0) = 1. 

The density contrast at the start of the matter-dominated era is determined by the primordial 
metric perturbation from inflation that we denote with £; n f. In Fourier space we have 


<Slin(k, Z) 


2 k 2 c 2 T(k)D(z) 


5SX 


* H o 


Cinf(k) . 


( 2 . 2 ) 


where the transfer function T(k) accounts for the damping of sub-Hubble-scale modes in the radiation 
era and T{k) —> 1 as k —> 0. Conventionally this is written in terms of a primordial Newtonian 
potential 

(5iin(k, z) = a(fc, 2 )d> in (k), (2.3) 

where, from eq. (2.3), we identify 




— - Cinf ) 
5 


v{k, z) = 


2 k 2 c 2 T(k)D(z) 


m, 


,H 2 0 


(2.4) 


The function a is plotted in fig. 1. 

Throughout this paper we will consider a primordial Newtonian potential with local-type non- 
Gaussianity, that is a local function of a Gaussian random field [71-73] 


•kin = <Pc(q) + /nl (^o(q) - (<Pg)) > 


(2.5) 


where v 5 G(q) is a Gaussian field, seeded by free field fluctuations during inflation, and /nl is a 
dimensionless parameter quantifying the magnitude of non-Gaussian corrections either due to non¬ 
linear evolution of the primordial metric perturbation [2], or to the non-linearity inherent in the 
general relativistic constraint equations [74-76]. Note that the non-Gaussian correction is a local 
function of the Gaussian field at a given initial (Lagrangian) position, q. 

In the presence of primordial non-Gaussianity, the linearly growing mode (2.3) contains first- 
and second-order terms with respect to tpc/q). It will be convenient, for our discussion, to define the 
first-order (Gaussian) density contrast 5c(k, z), where from eq. (2.3) and eq. (2.5) we identify 


(>G (k, z) = a(k, z)ipc (k). (2.6) 

The initial linearly growing mode of eq. (2.3) drives the subsequent non-linear evolution of the 
density field [79] and we may write the non-linearly evolved field in Lagrangian coordinates as 

5 = <5iin(q, Z) + <^onlin(q> z) . (2.7) 


where at second order we have 

<5no„lin(q,^) ^ ^n(q,*)) 2 + ^S 2 (q,z) (2.8) 

with s 2 = SijS l i and Sij, the trace-free tidal tensor , is defined as 

Sii= V- 2 £, (2.9) 

where dj ) is the Kronecker delta function. 
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Figure 1: The left panel shows a as a function of wavenumber k at redshift z = 0. The right panel 
shows the matter power spectrum P(k) at redshift z = 0. They are both obtained using CAMB [77] 
and assuming the Planck data set [78]. 


Since the Lagrangian and Eulerian coordinates (1.3) agree at leading order, the first-order density 
perturbation eq. (2.6) has the same form in either frame. However at second order the density 
perturbations differ because of the convective term proportional to the displacement, [80] 4 . For a 
general function f(x,z) we have, from eq. (1.3), 

/(x(q, z), z) ~ /(q, z) + 'S> ■ V/(q, z) . (2.10) 

Hence we can write the non-linearly evolved density field eq. (2.7) in Eulerian coordinates as 

S = <5i in (x, z) + C n ii n (x, z), (2.H) 

where at second order, using eqs. (2.8) and (2.10), we obtain [82]: 

Cnlin( X » z) ~ ~~ (<5lin ( x , z)f + ^S 2 (x, z) - ^(x, z) ■ VS(x, z) . (2.12) 

Similarly, although the Newtonian potential has the same form to first order, we must include a con¬ 
vective term when writing the primordial potential (2.5) up to second order in Euclidean coordinates 

v>g(x) ~ ^c(q) + ^(x, z) • v<£ G (q); (2.13) 

and thus 

$ in ~ ip G (x) - *( x , z) ■ V(pc( x ) + /nl (^g( x ) - (<Pg)) • ( 2 - 14 ) 

Hence we learn that, even though the primordial potential is by definition a local function of an 
initial Gaussian random field y>c(q) a t each initial spatial coordinate q and at a fixed initial time, the 
primordial potential at a fixed Eulerian position, x becomes time-dependent at second-order, due to 
the time-dependent relation between Eulerian and Lagrangian coordinates (1.3). In appendix A we 
show that even a Gaussian initial potential in the Lagrangian frame becomes a non-Gaussian field at 
second-order in the Eulerian frame due to the first-order displacement. We shall see that this gives rise 
to additional non-local terms in the distribution of collapsed halos in Eulerian space in the presence 
of primordial non-Gaussianity, assuming a bivariate local bias model (dependent on the linear density 
and the primordial potential) in Lagrangian space. 

4 See also [ 81 ] where this term is referred to as a shift term. 
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2.2 Peak-Background split 

Our local Lagrangian bias model relies on a simple though powerful argument, known as the peak- 
background split. In this approach, the first-order potential ip g is considered as a superposition of 
long and short modes 

¥>o(q) = y>G,z(q) + ¥>G,s(q), (2.15) 

which are statistically independent for a Gaussian random field. The local background is composed 
of long-wavelength modes ^>g,z and acts as an approximately homogeneous background cosmology on 
comoving scale l. On top of these the smaller scale peaks pg,s lead to the collapse of dark matter 
into halos, on a scale R -C l, when exceeding a suitable threshold value, S c . This is usually assumed 
to be the linearly growing density mode for a spherically collapsed object (S c = 1.686). 

With Gaussian initial conditions, this implies that the threshold for collapse is effectively different 
from place to place, 

S c —* S c - £> G ,z(q), (2-16) 

enhancing the formation of structures on top of long-mode overdense regions; this conclusion leads 
naturally to the idea that dark matter halos are biased tracers of the underlying matter distribution 
on a given scale l [83]. 

Substituting the peak-background split (eq. (2.15)) into the non-Gaussian primordial potential 
of eq. (2.5), we obtain 

$in(q) = <^G,Z + /nL (pg ,Z ~ (V^G,;)) + (1 + 2/nlV 3 G,z) <PG,s + /nL (Pg,s — Wg ,a]v) > (2-17) 

where the square brackets [p^ s ]y account for a local expectation evaluated over the volume V ~ / 3 
centred around q, 


[¥ J G,»]v(q) = 


lk>l 


dk 

-i (2^p 


!k‘>l 


dk' 
-i (2^)3 


3 -i(k+k').q^ G s(k)( ^ G s(k /)) . 


(2.18) 


By Fourier transforming eq. (2.17) and defining the Gaussian long and short density mode respec¬ 
tively as dim,G,z = cnpG,i an d <5ii n ,G,s = oip g , s ; we can identify a “background” density perturbation 

4n,z(k) = 8g,i + /nl« (pg,i — (‘Pgj)) ! (2-19) 

which changes the effective collapse threshold for the small scale peaks 

S c -^ 6 C - Vz • (2-20) 


At the same time, local-type non-Gaussianity introduces a correlation between long and short wave¬ 
length modes in the density field. This leads to the important conclusion that the small scale power 
is affected by the long-modes of the primordial potential, 


cq = (1 + 2 /nl^g, z)o"G • (2-21) 

The above discussion suggests that the number density of objects with mass M at redshift 2 in 
a volume V ~ I 3 , i.e. the local mass function (see appendix B for a brief introduction), depends not 
only on the mass and the redshift but also on the local density field and its moments [46, 84] 

n h = n h (M,z,[5? in ] v ). (2.22) 

In the non-Gaussian cosmology of eq. (2.5), eq. (2.22) for a given halo mass M and redshift z simply 
reduces to nh = ^h(^ii n> z, &i, /nl)- The dimensionless variable v in the mass function (eq. (B.4)) thus 
needs to be replaced with the local effective value 

v = ~ -► Kq) =- —rk -• (2.23) 

01 (q) 

This is the main result of this section and provides the basis for our local bias model in Lagrangian 
space. 
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3 Local Lagrangian bias model: halo overdensity 


Our bias model provides an expression for the local halo overdensity in Lagrangian coordinates 


$h( q) 


n h { q) - (n h ) 

(nh) 


(3.1) 


Using the preceding arguments we expect the local halo overdensity to be a function of the large-scale 
linearly growing mode of the density field, AirU (the local “background” density), and the small-scale 
variance of the linearly growing mode, cq (the “peaks”), averaged over the same large scale (l in the 
preceding section). We thus Taylor expand eq. (2.22) up to second order 3 in terms of Ain,; and cq 


- 1 


^(q) — AloAui ,1 + A)1 ( 

+ 2 Ao(^lin,z ) 2 + A) 2 (y— — 1^ + 2/3n<5lin ,1 ^ 


*-l 


(3.2) 


where we define the bias coefficients 


d l+J n h 

n h d^nn^ai 


>5lin,i=0,IT i —a 


(3.3) 


Hereafter we drop the subscript l, and simply use 5i; n and ipa to indicate the long-wavelength modes 
of the linearly growing density contrast and the Gaussian primordial potential respectively. 

The small-scale variance of eq. (2.21) depends on the local primordial potential (pc,i and hence 
we can write the Taylor expansion in the bivariate form [44, 67] 

£h(q) = ^ 10 Ain + froi^G + ^20(Ain) 2 + &uAin¥>G + &02^G . (3-4) 


(3.5) 


where we identify the Lagrangian bias coefficients 

&io = Aio , 

An = 2 /nlA)i , 

, L _ Ao 
&20 ~ 2 ’ 
b\i = 2/nlAli i 
^02 = 2 /nlA>2 ) 

Following BSS, we assume the Sheth-Tormen mass function corrected for non-Gaussian initial 
conditions (see appendix B for further details). This allows to provide explicit formulas for the 
Lagrangian bias coefficients; only two of them are independent: 

7^ 2 — 1 | 2 p 1 v 3 — v dn 3 /dM v + v~ x 

dln<r _1 /<iM 6 S c 


bio ~ 


pL __. a 7^-3 


b 20 = 'yv 


1 + ( 7 v 2 )p 8 C 
P 


k 3 - 


2 S r . 


+ 


(3.6) 


2<5? 


+ 


2 / yv 2 + 2p — 1 


1 + (7 p 2 )p 


S 2 


«3 

2 


'yp 5 - (7 + 2)p 3 + v 


81 


+ 


2 P 


1 + ('jp 2 )p S 2 


1 

+ - 


dn 3 / dM 


2 dlncr _1 /dM 


71/ 3 + (7 - l)is 

3 Jl 


+ 


2 p 


,-i 


1 + 3 8 2 


(3.7) 


where 7 and p were introduced in eq. (B. 6 ). Thus the remaining coefficients are obtained from the 
following combinations 

An = 2 /nlAA 10 > 

b\i = 2/ NL (A^2 0 — Ao), (3.8) 

bo 2 = 4 /nlA(AAo — 2 Ao) ■ 

In eq. (3.6) we can identify the first two terms as the usual Gaussian bias [85] plus a scale-independent 
correction introduced by PNG [10, 11], which has been showed to improve the comparison between 
theory and simulations [12]. The same structure is found for eq. (3.7). 

5 Note that to compute the tree-level bispectrum we need quantities up to 0(5 2 ). 
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4 Non-local Eulerian bias 


Galaxy surveys map the distribution of galaxies, which we assume to be located in collapsed halos, 
displaced with respect to their Lagrangian positions according to eq. (1.3). Equation (3.4) describes 
the excess of halos in Lagrangian space but this needs to be transformed to the Eulerian frame to 
account for their dynamics. 

Unlike the matter density contrast, the halo density contrast is not a 3-scalar since it is conven¬ 
tionally defined as a coordinate density [80]. Hence the number of halos in a given volume element is 
given by 

N h =n^(q,z)d 3 q = nf(x,z)d 3 x. (4.1) 


Thus we have 


1 + ( x , z) = [l + <S£(q,2)] 


d 3 q 
cPx 


(4.2) 


and using the coordinate Jacobian [see appendix D, eq. (D.5)] we obtain the transformation rule 

[65, 86] 

1 + (x, z) = [1 + <5(x, z)\ [l + <5£(q, ^)] , (4-3) 

The Lagrangian space halo density contrast, <5^(q, z), is given by eq. (3.4) in terms of the lin¬ 
early growing density contrast dun (q, z) and the Gaussian potential <pg(q) in Lagrangian coordinates. 
However, we have seen in eqs. (2.7) and (2.8) how <5ii n (q, z) can be expressed up to second-order in 
terms of the non-linear matter density, <5, and the tidal tensor, s 2 , 


4n(q, z) ~ <5(x, z) - ^(J(x, z)) 2 - =s 2 . 


(4.4) 


In eq. (2.13) we showed how ^(q) can be expressed up to second-order in terms of the displacement 
\l/ and the Gaussian potential in Eulerian coordinates </jg( x ) : 


^c(q) - </>g( x ) - 4f(x, z) ■ Vy> G (x) • 


(4.5) 


We thus obtain our expression for the Eulerian halo overdensity up to second order in terms of 
the density contrast and the Gaussian potential in Eulerian coordinates 

<^( x ) = bf 0 S + 6 0 V + 6 fo 52 + + bo2<PG - ^ b io s2 ~ & oi^( x > z ) ' V ^G ■ ( 4 -6) 

where we define the standard Eulerian bias coefficients (see fig. 2), 


b 

b 

b 

b 

b 


E — 1 -L 
10 — 1 ' 

E - h L 
01 ~ °01 


-'10 


20 


10 


= Tl b 
= bh + b\ 


e _ h L 
02 — ^02 ’ 


L 

20 


(4.7) 


Equation (4.6) generalises the result that is usually obtained under the spherical collapse approxima¬ 
tion (see appendix E). Indeed, the last two terms of eq. (4.6) are the non-local, non-linear terms of 
our bias model. While s 2 is an already known tidal term, we have derived here for the first time the 
convective contribution \t(x, 2 ) • V^q- We decide to keep the corresponding bias coefficients written 
as and frJb, instead of replacing them with bf 0 — 1 and 6 q X respectively. In this way we will be able 
to recognise more easily the differences they introduce with respect to the reference model of BSS, 
especially in the discussion in section 6. 

Finally, if we perform a Fourier transform with respect to Euclidean coordinates we obtain 

^( k ) = bf 0 6 + b^ipc + b% 0 6 *6 + bf ± S * </? G + b^VG * Vg - ^i 0 s 2 ~ 6 oi n2 . ( 4 -8) 
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Figure 2: The Eulerian bias coefficients as a function of mass and redshift, assuming /nl = 1- 


where we define 


dq 5 G (q)<5G(k-q) 


5 E (k) =S G {k) +f Nh a(k) f 

J {2nY a{q)a(\k — q|) 

« 2 ( k ) = / ^^5 2 (q,k- q)5 G (q)5 G (k- q) 


" 2(k)=2 /(l^ A/ ' 2(q,k_q) 

and the kernels are given by 


^G(q)fc(k - q) 
a(|k- q|) 


J ^l- 7r 2(q,k-q)5 G (q)(5 G (k-q) (4.9) 

(4.10) 

(4.11) 


^(ki, k 2 ) 

S 2 (k G k 2 ) 

A4(k G k 2 ) 


5 l k!-k 2 (h fc 2 \ 

7 2 k x k 2 \k 2 ki) 

(kr ■ k 2 ) 2 1 

k\k 2 3 
ki • k 2 
2k\ ' 


2 (ki ■ k 2 ) 2 
7 k\k 2 


(4.12) 

(4.13) 

(4.14) 


The standard second-order Newtonian kernel J- 2 is generated by the non-linear gravitational evolution. 
S 2 by the tidal term and the new kernel, A f 2 , is generated by the convective term M/ • Vy> G . 


- 9 - 


























isosceles 



k 3 /ki 


Figure 3: Explanation of the visual representation for the bispectrum introduced in [45]. The 
triangular-shaped region that hosts the colour map is due to the condition A : 3 < A ; 2 < k\. This 
requirement avoids double visualizations of the same triangular configuration. For the allowed values 
of k 2 /k\ and A: 3 /Aq we recognise same specific configurations: point (a) is for the squeezed limit 
(hi ~ k 2 A' 3 ), (b) for the equilateral configuration (Ay = k 2 = fc 3 ) and (c) for the folded one 
(k± = 2k 2 = 2fc 3 ). The elongated triangles (Aq = A: 2 + ks) resides on the left edge, while the upper and 
right edges correspond to isosceles triangles (k± > k 2 = A ; 3 or k-\ = k 2 > A; 3 ). General conhgurations 
are in the inner region. 


5 Three-point functions of halo and matter overdensities 

In sections 5.1 and 5.2 we will present the tree-level bispectra for the halo and matter overdensities. 
We adopt the definition 

(^(k 1 )<5|(k 2 )^(k 3 )) = (2Tr) 3 8 D (k 1 + k 2 + k 3 )S a ^(k 1 , k 2 , k 3 ), (5.1) 

where a, /3 ,7 = h, m with the labels h and m standing for halo and matter respectively. The crossed 
halo-matter bispectra are given in appendix G. 

In sections 5.2 and 6 we will make use of the graphical representation introduced by Jeong and 
Komatsu [45] to show the shape dependence of the bispectrum. The amplitude of B a / j 7 , or parts of 
it, will be plotted as a function of k 2 /k\ and k^/k\ in a colour map, under the condition < k 2 < Ay. 
This requirement avoids multiple visualizations of the same triangle. From all the possible choices of 
fc 2 /fci and jfe 3 /jfei, we can identify some specific configurations which are shown in fig. 3: equilateral 
(ki = k 2 = fc 3 ), isosceles (Ay > k 2 = k^ or ki = k 2 > fc 3 ), folded ( k\ = 2k 2 = 2k 3 ), squeezed 
(fci ~ fc 2 fc 3 ) and elongated (Ay = k 2 + A: 3 ). 

5.1 Matter bispectrum 

The matter bispectrum is given from eqs. (4.9) and (5.1) 

B mmm (k 1 ; k 2 ,k 3 ) = f2P(Aq)P(fc 2 ).T 2 (k 1 ,k 2 ) +2/ NL P(fcl \ P f 2 )Q; f 3) + 2 eye.) , (5.2) 

v a{k 1 )a(k2) ) 

where we have dropped the redshift z dependence from the function a(k, z) and the matter power 
spectrum P(k,z ) in order to simplify the notation; we will do the same in the following sections. 
The matter bispectrum is thus generated by both primordial non-Gaussianity, /nl, and non-linear 
gravitational evolution, P 2 . 

The tree-level approximation based on perturbation theory well describes the simulation results 
at scales up to A: ~ 0.05 — 0.1/iMpc - , depending on the redshift. Including one-loop corrections can 
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Figure 4: Shape dependence of the terms contributing to the halo bispectrum when Gaussian initial 
conditions are assumed and for k\ = 0.01,0.05,O.lft.Mpc^ 1 . Each term is normalized to the maximum 
value it can take in the A^/fcijfca/fci-space. This normalisation eliminates the redshift-dependence and 
as a result one should not make a comparison of the amplitude between plots because of the different 
scaling. Note that the last row shows the absolute value of M, since it can take negative values. The 
violet strip indicates where it is changing sign. 


significantly extend the validity for bispectrum to k ~ 0.3/iMpc -1 at redshift z > 1 [68]. Alternatively, 
it is possible to use an effective kernel .F 2 calibrated against simulations [87]. This phenomenological 
approach gives simpler expressions in the non-linear regime, and accurate predictions for the bispec¬ 
trum, up to k ~ 0.4/iMpc -1 in the redshift range 0 < z < 1.5, when Gaussian initial conditions are 
assumed [88]. 

5.2 Halo bispectrum 

First let us consider the halo bispectrum when Gaussian initial conditions (/nl = 0) are assumed [66] 

-®£ih( k i; k 2 > k 3 ) =b 3 w (2P(k 1 )P(k 2 )T 2 (k 1 ,k 2 ) + 2 cyc.) A 

+b 2 10 b 20 ( 2P(k 1 )P(k 2 ) + 2 cyc.) L (5.3) 

-^ b w b w (2P(fci)P(fc 2 )«S 2 (ki,k 2 ) + 2 cyc.) M . 

Here and in the following the various terms appearing in the halo bispectrum are labelled in accordance 
with BSS [44] in order to facilitate comparison and to highlight the new terms we have identified. The 
halo bispectrum of eq. (5.3) is generated by the non-linear gravitational evolution, P 2 , by non-linear 
bias, b 2 o, and, in particular, the last term, S 2 , is due to the non-local tidal term, s 2 , in the expression 
for the halo overdensity (see eq. (4.8)). Note that, here and in the following, we have suppressed the 
E superscript in the bias factors to simplify the notation but not the L superscript. This allows to 
keep track of the effects generated by the non-local terms s 2 and n 2 . 

In fig. 4 we plot the shape dependence of each of the terms A,L,M for different choices of fci, 
normalizing each of them to the maximum value it takes in the {k 2 /ki, fc 3 /fci)-space. This normali¬ 
sation eliminates the redshift-dependence and as a result one should not make a comparison of the 
amplitude between plots because of the different scaling. The results are shown in fig. 4. Note that 
the absolute value is plotted for M, as it can be positive or negative. Indeed, the violet region that 
cuts the plots into two parts is where M changes sign. 
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Figure 5: The term A is given by the sum of two pieces, which we label Ai and A 2 . The former 
is sourced by non-linear gravitational evolution while the latter by PNG. We show them separately, 
normalized to the maximum value A can take and assuming z = 0 and /nl = 10. The normalization 
does not completely cancel the redshift dependence, which is present in A 2 through the factor D(z)~ x . 
The effect of A 2 is visible in the squeezed configuration, where Ai is vanishing and A 2 is at its 
maximum (the logarithmic scale might hide this aspect but a quick look back to fig. 4 should clarify 
this point). 


In the presence of primordial non-Gaussianity, many more terms contributes to the halo bispec¬ 
trum: 


-Bhhh(ki, k 2 , k 3 ) — ^hhh >L '(k 1 , k 2 , k 3 ) 


~jb w b 10 (2P(fci)P(fe 2 )5 2 (ki,k 2 ) + 2 cyc.) M 

-^ 10601^0 (mki)p(k 2 ) + ~^ jr 

7 \ VQWi) a ( k 2) 


<S 2 (ki,k 2 ) + 2 eye. 




2 iL 


7 


01 w 10 


a{ki)a{k 2 ) 


-bioboi ( 2P(fci)P(fc 2 ) 


A/" 2 (ki,k 2 ) A/2(k 2 ,ki) 


-& 10 & 01&01 2P(fci)P(fc 2 ) 


a(k 2 ) cx(ki) 

A7 2 (k!,k 2 ) , M(k 2 ,ki) 


+ 2 eye. 


a(k 2 


ot{k 1 ) 


1 




^(fcQPfe) ( N 2 {^M) + A/~ 2 (k 2 , k 3 ) + 2 cyc 

ot(ki) 


a(ki)a(k 2 ) \ a{k 2 ) 


1 


(5.4) 


ot{ki) a{k 2 ) 


+ 2 cyc. 


The quantity Pi 1 hh >L '* accounts for the terms with label going from A to L; it matches exactly Eq.(5.1) 
of BSS, that we reproduce in appendix F. 

The first line in eq. (F.l), term A, comes from linear bias acting on the matter bispectrum. We 
can split this into two terms 


Ai(ki,k 2 ,k 3 ) =2P(fci)P(fc 2 )P 2 (ki,k 2 ) + 2 cyc. 


A 2 (fci, fc 2 , fc 3 ) =2/nl 


P(fci)P(fe 2 )a(fc 3 ) 

a(k\)a{k 2 ) 


+ 2 cyc.. 


(5.5) 

(5.6) 


identifying an additional term, A 2 , with respect to the Gaussian initial conditions eq. (5.3) propor¬ 
tional to /nl- We study the shape dependence of A 1 and A 2 in fig. 5. Again, we plot them separately 
but normalize to the maximum value taken by A(ki,k 2 ,k 3 ) = Ai + A 2 in the A: 2 //ci,fc 3 /fci-space. 
The relative values of the two plots can thus be compared, but note that the normalization does not 
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Figure 6: Shape dependence of the terms, B to F, contributing to the halo bispectrum. These are 
generated by non-Gaussian initial conditions. A value /nl = 10 and redshift 2 = 0 are assumed. 
Each term is normalized to the maximum it can take. This choice does not completely cancel the 
redshift dependence in the terms B and C, where it is present as a factor D(z)~ 1 , and does not 
allow a comparison of the amplitude between different plots. We clearly see that the effect of PNG is 
prominent in the squeezed configuration. 


cancel the growth factor 1/a oc 1 /D(z) in A 2 but absent in Ai. Hence the relative amplitude of 
the bispectrum from primordial non-Gaussianity A 2 grows with redshift relative the bispectrum Aj 
coming from Gaussian initial conditions. Figure 5 highlights the interesting shape dependence of A 2 ; 
it peaks in the extremely squeezed configuration (top left), exactly where the Ai term vanishes (see 
also fig. 4). 

We emphasize that A 2 and the terms from B to K are generated by primordial non-Gaussianity. 
We plot the shape dependence of the terms B to F in fig. 6 and the term G to K in fig. 7. They 
confirm what we have previously stated: the PNG terms are greatest in the squeezed configuration. 

As a result of the presence of s 2 and our new term n 2 in the expression for <5® (see eq. (4.8)), 
additional terms appear with respect to BSS. M, N and O that are generated by the tidal term s 2 ; 
schematically, they are sourced by ( s 2 S5 ), ( s 2 5tp ) and (s 2 ipip) respectively. We have seen in eq. (5.3) 
that M is present regardless of the presence of PNG, but N and O come from the coupling between 
the tidal term and ipQ for / NL ^ 0. On the other hand, P, Q and R are generated by n 2 . They 
account for the non-local effect of the potential ipQ and are therefore due to the presence of PNG. 
Schematically, they are generated by ( n 2 55 ), ( n 2 6<p ), ( n 2 iptp ) respectively. In fig. 8 we show the terms 
from N to R. Since they can take negative values, we plot their absolute value. Interestingly, in all of 
the plots we can identify a violet strip, indicating a change of sign and, hence, where the kernels 5 2 
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Figure 7 : The terms contributing to the halo bispectrum, with label going from G to K. These are 
generated by non-Gaussian initial conditions. A value /nl = 10 is assumed. Each term is normalized 
to the maximum it can take. This choice completely cancels the redshift dependence and does not 
allow a comparison of the amplitude between different terms. We clearly see that the effect of PNG 
is prominent in the squeezed configuration. 


and A4 make each term vanish. 


6 Analytic estimates 


In this section we further investigate our result for the halo bispectrum and, in particular, we compare 
it to our reference model BSS [44]. Since the model of BSS has shown a good fit against simulations 
[44, 89], we want to understand if and where differences between these two models arise. 

In fig. 9 we plot the absolute value of the relative difference between our halo bispectrum and 


that of BSS 


Diff(k|, k 2 , k 3 ) 


l?TRTw(ki, k 2 , k 3 ) — -E>BSs(ki, k 2 , k 3 ) 
-Bbss^i, k 2 , k 3 ) 


( 6 . 1 ) 


for values of k\ = 0.01,0.05, CI.lft.Mpc~ 1 and redshift z = 0,0.5,1. We consider halos of mass M = 
10 13 1i -1 Mq and primordial non-Gaussianity with /nl = 10. We choose to saturate differences above 
10 % to the same red colour of the palette for the purpose of presenting many different plots with 
different ranges of values in a compact way. However we do not observe differences above 25% 6 . 

6 An exception is the top left plot of fig. 9, where there is a small curve close to the squeezed configuration in which 
the discrepancy is actually bigger than that. This because in that area the BSS halo bispectrum crosses zero for these 
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Figure 8: Shape dependence of the terms with label going from N to R contributing to the halo 
bispectrum. These are the additional terms generated by the presence of s 2 and n 2 in 5h- A value 
/nl = 10 is assumed. Each term is normalized to maximum value it can take so that the redshift 
dependence is completely dropped; note however that the different scaling does not allow one to 
compare the amplitude between plots. Since these terms can take negative values, we show their 
absolute value. Actually, the blue line indicates where they change sign, with the top right part of 
the plots being positive. 


For the value of the halo mass and /nl considered, we find that the terms M and P are the 
main sources of the differences, with all the other terms contributing very little or having a negligible 
effect. The most relevant differences appear for k\ = O.Ol/iMpc^ 1 : up to 25% in the elongated, folded 
and equilateral regions for all the redshifts considered. These discrepancies drop to a few percent for 
k± = 0.05/iMpc^ 1 at z = 0, while for z = 0.5,1 they are reduced to about 5% when approaching the 
equilateral configuration and to order 10% in the elongated and folded regions. For k\ = 0.1/iMpc -1 
we observe a similar pattern, but at z = 0.5 the approximately 10% difference area that was present 
in the elongated and folded regions for k\ = 0.05^Mpc _1 is almost completely washed out, decreased 
to about 5%. However, that area is still present for z = 1, although reduced in size. We also recognise 
an area where the difference is close to zero [the violet region going approximately from the squeezed 
(top left) to the isosceles (bottom right) configuration] corresponding to a vanishing contribution from 
the terms M to R, as shown in figs. 4 and 8. Interestingly, in all plots the squeezed configuration is 
unaffected. 

We can understand the features that we have just described by studying the analytic solutions 
of the halo bispectrum in three simple configurations: 

values of the bias coefficients. 
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Figure 9: The relative difference in absolute value between our model and the one by Baldauf et ah, 
assuming M = 10 13 /i _1 Mq and /nl = 10. Note that we saturate differences above 10% to the same 
red colour of the palette; as the plots range between different values, this choice allows us to present 
them in a compact way, highlighting where the most relevant differences are expected. However we do 
not observe discrepancies above 25%. Interestingly, the squeezed limit is not affected while the other 
configurations show differences from only a few percent. A no-difference region going approximately 
from the squeezed to the isosceles configurations is present, following the shape dependence that we 
showed in figs. 4 and 8. 


In the equilateral configuration, where k\ = Aq = Aq = k, the halo bispectrum becomes 


b 2 

Shhh =\ vjr(126lO + &10 + 42620 ) 


1 

a{k) 


7 boibio(l2bio + b\ 0 + 426 2 q) + 36f 0 ( 6 q X + 2(/nl6io + ^ 11 ) 


1 

\ b2 o 1 , 

a(k) 2 

7 

3 

L 2 [ 

a(k) 3 

°01 1 


(l 26 io + &10 + 42 & 2 o) + 6601610 ( bg X + 2 (/nl6io + 6n) ) + 66026 

6&01&02 


&01 ( &01 + 2(/ nl&io + bn ) ) + 84601602610 


+ 


x(k) 4 


>p{ky 


( 6 - 2 ) 


We remind the reader that the contributions of b\ 0 and 6^ indicate where the non-local, non¬ 
linear terms of our model (s 2 and n 2 ) are introducing differences respect to the BSS model. By 
looking at eq. (5.4) we notice that the terms M, N, O, P, Q, R are respectively linked to the bias 
combinations 6 2 0 6(' 0 , 6io6oi6^ 0 , ^oi^ico ^io^oii ^ 10 ^ 01^017 ^oi^oi> so that we can actually recognize 
each of them in eq. (6.2); M appears in the first line of eq. (6.2), while N and P appear in the 
second line, O and Q in the third and R in the first term on the fourth line. 

The presence of a larger difference at high redsliift and on large scales can be explained by noting 
that the function a(k) takes smaller values in those regimes and enhances the terms inside the 
square brackets (see eq. (2.3) and fig. 1) and, therefore, accentuates the effect of b\ 0 and 6^. 


• In the folded configuration, &2 = k% = k and Aq = 2fc, the halo bispectrum can be written as 


I?hhh — 


26? 0 (2 — n 2 ) + 6 2 0 



(1 + 4n 2 ) + 
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where we define II = T(2k)/T(k) to make the notation more compact. The full expression is 
long so we have given just the first two terms. These are enough to explain what we observe 
in fig. 9. Again, b\ 0 and 6^ are present. As for the equilateral configuration, the effect of M 
appears in the first term, while N and P in the second one, O and Q in the third and R in the 
fourth one. We see that the effect of b\ 0 and bg 1 can be enhanced depending on the value of 
II. When k\ = 0.01/iMpc the ratio II « 1, while for the other two values of k\ the ratio II < 1. 
This suggests why bigger differences should be expected in the case with k\ = 0.01/iMpc. 

The same considerations apply to the function a(k ) as in the equilateral case, so that the 
differences are larger at high redshift and on large scales. 


• A simple expression for the halo bispectrum in the squeezed limit can be found by setting 
k\ = &2 = ekg = k with e> 1, corresponding to an isosceles triangle, whose degree of squeezing 
is controlled by the parameter e. For squeezed triangles (large values of e), the leading term in 
the bispectrum is 


/ii 


hhh 


a(k)- \ 


( 2/nL&01^10 + bgibigb 


-'ll 


^ W b2 ° lb ° 2 


P(k) 


2 .3 


l(fc ) 3 


2/ nl ^ oi^io + &0A1 + 2601^02^10 


(6.4) 


The absence of terms involving bias coefficients b\ 0 and bg 1 explains why our predictions do not 
differ from the BSS model in extremely squeezed configurations. 


7 Conclusions 

An important application of measurements of large-scale structure in our Universe is to determine the 
distribution of primordial perturbations, in particular possible non-Gaussian signatures of scenarios 
for the origin of structure in the very early universe. For example, the shape information contained 
in the primordial bispectrum is a valuable tool to discriminate between different inflationary models. 
The matter bispectrum at later times is due to a combination of both primordial non-Gaussianity and 
non-linear evolution under gravity. However, the way in which gravitationally collapsed halos trace 
the density field, the bias model, can enhance the effect of primordial non-Gaussianity in the galaxy 
distribution. In particular, local-type non-Gaussianity leads to a scale-dependent bias which can have 
a dramatic effect on very large scales in both the halo power spectrum [8, 9] and the halo bispectrum 
[!«>:• 

In this paper we developed a local Lagrangian bias model, focussing on second-order, non-local 
and non-Gaussian effects. In particular, we extended and applied a local Lagrangian biasing scheme 
to a general set-up with local-type primordial non-Gaussianity. For an /nl cosmology, we re-derive 
the known result [67] that the halo overdensity can be expressed as a bivariate expansion in terms of 
the linear matter overdensity in the Lagrangian frame and the primordial Gaussian potential. 

Non-linear evolution of the matter field in general gives rise to second-order terms in the matter 
density which include non-local, tidal terms (see eq. (2.8)), while transforming from the Lagrangian to 
the Eulerian frame introduces a non-local convective term at second order (see eq. (2.12)). Non-local 
here implies terms derived from derivatives of the potential, not directly from the local density or 
its derivatives. Both these terms are included in the usual kernel, Ti, for the second-order density 
in Eulerian space. But since the halo density is determined by the linear matter overdensity in the 
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Lagrangian frame, one must account separately for these non-local terms to reconstruct the halo 
density at late times. These terms are absent at early times, or if we restrict ourselves to the spherical 
collapse approximation (see appendix E). 

We have shown in this paper that in the bivariate expansion we must also account at second 
order for the convective term relating the primordial potential in the Lagrangian frame to that in 
Eulerian space at later times (see eq. (2.13)). This gives rise to a new term in the halo bispectrum 
in the presence of local-type primordial non-Gaussianity which has not previously been studied as far 
as we are aware. 

Setting / NL = 0 and are able to recover the halo bispectrum model of [ 66 ], when a local La¬ 
grangian biasing scheme is applied. Three terms appear in the halo bispectrum (A,L and M) that are 
sourced by (A) the non-linear matter density encoded in the kernel J~ 2 , (L) the non-linear bias 620 
and (M) the tidal term s 2 . 

Generalising to /nl 7^ 0, we found 12 terms in the halo bispectrum (labelled A to L) matching 
the BSS model [44]. In this case, the non-linear matter density term, A, includes a correction due 
to PNG, while the non-linear bias term, L, is left unchanged. The other contributions (B to K) 
come from a mixture of bivariate terms, involving the bias coefficients &oi 5 & 11 , t> 02 - We also found a 
contribution, M, sourced by the tidal term, s 2 , which also couples with terms in the halo overdensity 
that are specifically due to PNG and, hence, generate new contributions, N and O, in the halo 
bispectrum. The new convective term, n 2 , also generates contributions in the presence of primordial 
non-Gaussianity; in the halo bispectrum we have found three new contributions, P, Q and R, due to 
this term. 

In order to investigate the magnitude and shape of the various contributions to the bispectrum, 
we have implemented a version of the Sheth-Tormen mass function corrected for PNG in light of 
the Lo Verde mass function, following BSS [44]. This allowed us to numerically calculate the bias 
coefficients for our model and predict the halo bispectrum in different configurations for sample values 
of /nl and at various scales and redshifts. 

We investigated the halo bispectrum by comparing it to the fiducial model of BSS. Assuming 
halos of mass M = 10 13 1 i _1 Mq and /nl = 10, we found: 

- At redshift z = 0 differences up to 25% in the halo bispectrum for k\ = 0.01/iMpc -1 in the 
elongated and folded configurations and approximately 7 — 8 % when approaching the equilateral 
configuration, while these drop to a few percent when k\ = 0.05 or O.lIiMpc^ 1 . 

- At redshift z = 0.5 differences up to 25% for k\ = 0.01/iMpc -1 in the elongated, folded and 
equilateral configurations. When £q = 0.05/iMpc _1 differences of order 10% are still visible in the 
elongated and folded shapes, while they decrease to approximately 5% towards the equilateral 
configuration. For fci = 0.1/iMpc -1 these differences reduce to about 5%. 

- At redshift 0 = 1 we find results similar to those for 0 = 0.5, except that differences of order 
10% are still visible in the elongated and folded regions for fci = 0.1/iMpc -1 . 

In general, we observe that the non-local terms have a negligible effect in extremely squeezed config¬ 
urations. 

Our results indicate that the non-local terms in the halo overdensity could have a significant 
contribution to the galaxy bispectrum, especially on large scales and at high redshift. The next 
challenge is to test these theoretical predictions against N-body simulations with non-Gaussian initial 
conditions. Ultimately we would wish to be able to estimate the signal-to-noise for upcoming surveys, 
like the ESA Euclid mission [90], which probes large scales and high redshifts, in order to explore 
the observability of these non-local effects in the bispectrum. A full discussion of the observability of 
these effects must include a halo occupation model, to describe how galaxies populate halos and many 
other effects including redshift space distortions and lensing along the line of sight, from galaxies to 
the observer, in order to translate theoretical model for the halo bispectrum into predictions for the 
observed galaxy angular bispectrum in redshift space. Nonetheless, we have identified in this paper 
novel contributions to the expected galaxy bispectrum on large scales and high redshift in the presence 
of primordial non-Gaussianity. 
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A Displaced Gaussian random fields 

In deriving the halo abundance in Eulerian space we need to map the initial gravitational poten¬ 
tial in Lagrangian coordinates, eq. (2.4), into Eulerian coordinates under the coordinate displace¬ 
ment eq. (1.3), which is itself determined by the gravitational potential. In this appendix we shall 
demonstrate how even an initial Gaussian held in Lagrangian coordinates may be transformed into a 
non-Gaussian held by this displacement to Eulerian coordinates. 

We start from a random held, <^>g( ( ?)i defined with respect to the Lagrangian coordinate chart, 
q. Let 

<PG(q) = e/(g)d, (A.l) 

where a denotes a Gaussian random variable and e is a small perturbative parameter. Thus <Pg(<i) is 
a Gaussian random held (hrst order with respect to e) with, for example, vanishing 3-point function 

{<PG(qi)0G(q2)0G(qa}) = e 3 /(gi)/(g 2 )/(d 3 )(a 3 ) = 0 . (A.2) 

where angle-brackets denote the ensemble average. 

Let the Eulerian coordinate x be related to q by a hrst-order displacement held ^(g), correlated 
with the held (pc such that 

x(q) = q +e%/j{q)a. (A.3) 

If we consider a hxed coordinate q then x(q) is itself a random variable, correlated with pg{q)- We 
can then construct the held 


<Pa(x{q)) = <p G (q) + e 2 ip(q)f'(q)a 2 + 0(e 3 ), (A.4) 

which is Gaussian at hrst order in e, but non-Gaussian at second order. For example, the 3-point 
function of <p(x(q)) with respect to the coordinate chart q is non-vanishing at fourth order 

{<PG(x(qi))<PG(x(q2))<PG(x(q 3 ))) q = e 4 [-tf’(qi)f'(qi)f(q 2 )f(q3) + perms] (a 4 ) + 0(e 6 ) ■£ 0 

Conversely, if we work with respect to the Eulerian coordinate chart x, it is the Lagrangian 
coordinate that becomes a random held at hxed coordinate x: 

q(x) = x — eijj(x)a + 0(e 2 ), (A.5) 

and Lpc{q{x)) becomes a non-Gaussian held at second order 

l PG(q(x)) = <Pg(x) - e 2 ip(x)f(x)a 2 + 0(e 3 ) , (A.6) 

where <Pg(x) = ef(x)a is a hrst-order Gaussian random held in Eulerian space. 

B Halo mass function 

The Press-Schechter formalism [91] and its extensions [92, 93] provide a model to describe the full non- 
linearly evolved density held, and in particular the number of gravitationally collapsed dark matter 
halos, in terms of the initial, linearly growing density held (see [94] for a pedagogical review). Dark 
matter halos are identihed as peaks in the linearly growing density held of eq. (2.3), exceeding a 
suitable threshold value, S c . This is usually assumed to be the linearly growing density mode for 
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a spherically collapsed object. For example, in a flat ACDM universe the threshold for a spherical 
collapsed halo is [95] 

S c (z) = 3 ( 2 2 q 2/3 [1 + 0-0123log , (B.l) 

reducing to S c ~ 1.686 when S7 m = 1. Hence S c is weakly dependent on the value of Q m and Ha at 
the time of collapse [96]. 

For Gaussian initial conditions, the smoothed first-order density field is a Gaussian field with 
variance 

*g(M,*0 = / dkk 2 W^(k,R)P 0 (k), (B.2) 

where Po{k) is the linear matter power spectrum at redshift z = 0 (see fig. 1) and Wm(/e, R) is a window 
function in Fourier space that filters out modes below the length scale R(M) = (3M/ 47 r Pm) 1 ^■ We 
will adopt the real-space top-hat filter, with Fourier transform 


W M (k, R) 


3 

(Ml ) 3 


[sin(fcR) — fcRcos(fcR)] . 


(B.3) 


The Press-Schechter (PS) approach predicts the number density of objects with mass M at 
redshift z, i.e. the mass function , to be 


n h (M,z) 


/M 


Pm 

"m 


din <7 g 
dM 


(B.4) 


where we introduce the variable v = 8 c /aQ. 


The analytic form for the PS distribution is 


/psM = 



(B.5) 


By considering the collapse of ellipsoidal overdense regions, the Sheth-Tormen (ST) distribution 
[97-99] is obtained using 


/st(p) = Mp) 



~P 


ve 


-7 


u 2 

~2~ 


(B. 6 ) 


where the parameters 7 = 0.707 and p = 0.3 are found by fitting against simulations and A(p) = 0.322 
under the requirement that all the mass is collapsed into halos. Equation (B. 6 ) greatly improved the 
agreement with simulations. Interestingly, both the PS and ST functions, /, depend only on the 
variable v, for this reason they are know as universal mass functions. 

In presence of weakly non-Gaussian initial conditions, the probability distribution function (PDF) 
of fluctuations can be approximated by an Edgeworth expansion. Then, a derivation similar to the 
one yielding to eq. (B.5) leads to [100] 


/lv(g M) = /ps(p) 


1 + ^ ^k 3 (M )H 3 (v) 


dK 3 (M)/dM H 2 (v)\ 
dlncr -1 /dM v ) 


(B.7) 


H n is the n-th Hermite polynomial and k 3 (M) the 3rd cumulant, defined as k 3 (M) = (^)/a 3 , where 

( S L) = J J ( 7 ^ J ^^W M {p)a{p,z)W M {p')a{p',z)W M {p'')a{p\z){^ in {p)^ in (p>')^- m {p,'')). 

(B- 8) 

In addition, for primordial non-Gaussianity of the form of eq. (2.5), the variance needs to be replaced 
with 


= < 4 „>=/ J -^^WM(p)a{p,z)W M {p , Mp , ,z)(^ in {p)<S>i n {p 1 )} 


(27r) 

CTq (1 + «2(M)) 


(B.9) 
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Figure 10: The left panel shows the variance cr at redshift z = 0. Note that values of /nl inside the 
current constraints from CMB and LSS produce no appreciable effects on a. On the right side the 
2nd and 3rd cumulant and the derivative of the latter are shown, assuming /nl = 1. 


Equation (B.7) is known as the Lo Verde et al. (LV) mass function. In [101] fitting functions for k 2 
and K 3 are given; although k 2 oc /^ l , it gives a negligible correction to cr a for any realistic value of 
/nl and we will neglect it here 7 . The 3rd cumulant K 3 reads 


k 3 (M) 


/nl ( 6.6 x IQ" 4 ) 


1 - 0.016 In 


M 


h - 1 M, 


0 /J 


(B.10) 


we refer the reader to appendix C for further details on K 3 and 

Note that the LV mass function is no longer universal, although the explicit dependence on the 
mass M through eqs. (B.9) and (B.10) is weak (see fig. 10). Throughout the paper we follow BSS and 
use an effective form of the LV mass function by replacing fj,v with 

/lv —» /lv 7^- ; (B.ll) 

J PS 

the resulting mass function can be thought as the ST mass function corrected for non-Gaussian initial 
conditions. 


C Comment on the fitting functions 

In [103], useful fitting functions for K 3 and d ^*1 are given 

4 1} = 0.000329 (1 + 0.09z) &/ 0 09 

d^/dM = _ 000Q61 0 . 25 

dhicr-i/dM v ' 1 

where b± = (v 2 — 1 )/ 6 c + 1 is the Eulerian Gaussian bias derived from the PS mass function with 
5 C = 1.42. 

'In general re 2 °c r [(-i J //\| j - The model of eq. (2.5) satisfies the Suyama-Yamaguchi equality: T\| = (6/5) 2 /^ L [102], 
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As we use a combination of the ST and LV mass function with S c = 1.686, we prefer to use 
eq. (B.10) for K 3 and get out of this one. We obtain 


d,K3/dM 

dlner -1 /dM 


1.056 x 10 -5 /nl 


4.2xlO~ 10 j-2 
1 + K 2 ■'NL 


/nl 1 


(C.l) 


where we introduced the quantity Tq that is defined as 


dk 

(2tt) s 


W M jo{kR)P 0 (k ), 


(C.2) 


and jo(kR) is the spherical Bessel function of order 0. K 3 and g 3 -^/dM are plotted in fig- 10. 


D Coordinate Jacobian 

Since the matter density is a 3-scalar, we have 

M = p(x, z)a(z ) 3 y/\g^\d, 3 x = p{ q, z)a(z ) 3 y/\g^\d 3 q , 


(D.l) 


the determinant of the Eulerian metric g x is simply \g x \ = 1, but the Lagrangian space has a non-trivial 
metric g q even in Newtonian theory. If we define the coordinate Jacobian 


J = 


d 3 3 


d 3 q 


= V 1 9q 


by using eq. (D.l) we have to first order 


3 ~ 1 + V-*. 


(D.2) 


(D.3) 


At the initial redshift Z; n , the Eulerian and Lagrangian frame are equivalent, i.e. \I/ = 0, and hence 
J = 1. If in addition we assume that the mass per volume element is conserved and the initial density 
was uniform, p(q, 2 ; n ) = p{z i n ) in the limit z- ln 00 , we have 


M = a 3 (z)p(x , z)d 3 yi = a 3 (z in )p(z in )d 3 q 


and hence 


J = 


d 3 ) 


d 3 q 


a 3 (z in )p(z in ) p(z) 


a 3 {z)p{yi 1 z) p(x,z) 


= [1 + J e (x, 2 )] 


-1 


(D.4) 

(D.5) 


where J E (x) is the Eulerian fully non-linear density contrast. 


E Spherical collapse approximation: Local Eulerian biasing 

The spherical collapse approximation has been used in many works (see for instance [44, 63, 67]) and 
we present it here with the aim of highlighting the differences respect to the more general result we 
derived in section 4. 

The spherical collapse corresponds to the special case in which tk = 0 and the velocity field 
vanishes at the centre of the symmetrical collapse so that the Eulerian and Lagrangian coordinates 
coincide at all times, x = q. The results that we are going to quote all relate the density at the 
centre of collapse or assume a uniform density equal to that at the centre, i.e. edge effects are not 
considered. 

Following [63] (see also [65]), the linearly growing density at the centre of the collapse relates to 
the non-linear density field 


^Un(q) = + «2^ 2 + ■ ■ - , (E.l) 

i=o 
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where the coefficients are 


ai = 1 , a2 = ~2i' ( E - 2 ) 

Equation (E.l) is local in either Eulerian or Lagrangian coordinates since they coincide at the centre 
of the symmetrical collapse. 

On the other hand, at the centre of collapse the auxiliary potential is simply 

¥>G(q) = ¥>g(x), (E.3) 

where we remember that <^g( x ) always refer to the primordial value as it simply allows to keep track 
of the PNG effects in the initial density field; hence </?g( x ) does not evolve with time like a physical 
quantity do and its coordinate transformation is straightforward. 

Now, using eqs. (E.l) and (E.3) to express <5£(q) = <5^(q(x, r)), we can solve eq. (4.3) for <$®(x). 
By using the same bias redefinition of eq. (4.7), the Eulerian halo overdensity finally reads 

S h (x) = bf 0 d + b^ip G + b% 0 5 2 + b^Sipd + b^ 2 G . (E.4) 

We see that the local Lagrangian biasing of eq. (3.4) is compatible with the local Eulerian model we 
have just obtained, despite the fact that the transformation of eq. (4.3) is inherently non local: this 
is true only when spherical collapse dynamics applies. 


F BSS halo bispectrum 

The halo bispectrum model of BSS [44] predicts 


kr,k 2 ,k 3 ) =bl o (^2P(fc 1 )P(fc 2 )J r 2 (k 1 ,k 2 ) + 2/ NL P M^1 + 2 eye. 

V a{k\)a(k 2 ) 


+b 2 10 b 01 2 P(fci)P(fc 2 ) 


2/nl 


1 


1 


a(fci) a(k 2 ) 

1 1 


F 2 {ki,k 2 ) 


+ 2 eye. 


&oi 


+ ^10^1 
+ ^ 01^02 ( 2 


P(fci)P(fc 2 )q(fc 3 ) 

a(ki)a(k 2 ) \a(ki) ot{k 2 ) 

P(fci)P(fc 2 ) , ^ , OJ , P(fci)P(/c 2 )a(fc 3 ) 

2 . —Tw^p2(ki,k 2 ) + 2 / N l— _ 9/ ,_ x _ 9/ .,_ ^ -1- 2 eye. 


a(ki)a(k 2 ) 

F(fci)P(fc 2 ) 


+boib w b 02 2 


a 2 {ki)a 2 {k 2 ) 
P(ki)P(k 2 ) 


+ 2 eye. 
1 


a(ki)a(k 2 ) \a(ki) a(k 2 ) 


a 2 (ki)a 2 (k 2 ) 


+ 2 eye. 


+^10^02 2 


P(fci)P(fc 2 ) 


+ 2 eye. 


a(ki)a(k 2 ) 

+b 2 01 bn ( P S i Y7FT f ■ ~lk) + + 2 cyc - 

\a(ki)a(k 2 ) \a(ki) a(k 2 )) 


(F.l) 


+boibi 0 bn ^P(fei)P(fe 2 ) 

+b 2 w bn (p{k\)P{h 2 ) 


1 


1 


a 2 (ki) a(k\)a(k 2 ) a 2 (k 2 ) 

1 1 


+ 2 cyc. 


+^01^20 ( 2 


,P(fci)P(fc 2 ) 


a(ki)a(k 2 ) 
+boibiob 2 o f 2P(ki)P(k 2 ) 


a(k i) a(k 2 ) 

+ 2 cyc. 

1 


+ 2 cyc. 


1 


a(k i) a{k 2 ) 


+ 2 cyc. 


+b 2 10 b 20 (2P(ki)P(k 2 ) + 2 cyc.) L . 


Equation (F.l) is incorporated in our prediction for the halo bispectrum but new terms appear (see 
eq. (5.4)). 
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G Halo-matter bispectra 


Here we consider the case of crossed bispectra between halos and matter. Potentially, weak lensing 
measurements will allow to cross correlate the dark matter density field with galaxies in the future. 
Also, these results provide additional predictions of our model that can be tested against simulations. 
We start by quoting the halo-halo-matter bispectrum in presence of Gaussian initial conditions 


-Bhhm(ki,k 2 ,k 3 ) =b\ o 6J‘2(k 1 ,k 2 )P(fci)P(fc 2 ) + 2 eye. 


+610620 4 P{k 1 )P(k 2 ) + 2 eye. 


(G.l) 


-^6i 0 6i 0 ( 25 2 (k 1 , k 2 )P(fci)P(fc 2 ) +2 eye. 


Again, we recognise that it is sourced by non-linear gravitational evolution and non-linear bias, while 
the last term is generated by s 2 in eq. ( 4 . 8 ). 

Assuming PNG, the halo-halo-matter bispectrum 8 reads 


Phhm(ki,k 2 ,k 3 ) =6 2 0 ( 6P 2 (kx, k 2 )P(fc!)P(fc 2 ) + 6/ NL a(fc 3 ) + 2 eye. 


+6 01 6 10 f4P 2 (ki, k 2 )P(fc 1 )P(fc 2 ) 
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(G.2) 


+ 2 eye. 


The terms with label going from A to I match exactly the result of Eq.( 5 . 6 ) in BSS but, as for the halo 
bispectrum, new terms appear. J, K are due to the tidal term s 2 : schematically, they generated by 

®Note that the term H corrects a typo present in Eq.(5.6) of BSS. 
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(s 2 <5W<5«), (s 2 6 ( ' 1 ' > ip) respectively. We see by comparison with eq. (G.l) that J is present regardless of 
PNG, while K comes from the coupling between the tidal term and ip which is specifically introduced 
by PNG. L and M are present because of n 2 and, therefore, depend on the presence of PNG. These 
are schematically generated by (n 2 ^ 1 )^ 1 )), ( n 2 5^ip ), respectively. 

Finally, for the halo-matter-matter bispectrum when Gaussian initial conditions are assumed we 

find 


£hmm(ki,k 2 ,k 3 ) =&io 6J r 2 (ki,k 2 )P(fci)P(fc 2 )+ 2 eye. 


+b 20 ^2P(fci)P(fc 2 ) + 2 eye. 
-y&i 0 ^25 2 (ki,k 2 )P(fci)P(fc 2 ) + 2 eye. 


(G.3) 


where the same considerations for eq. (G.2) apply. Then, assuming PNG, the halo-matter-matter 
bispectrum reads 

Phrara(ki,k 2 ,k 3 ) =& 10 f6P 2 (k 1 ,k 2 )P(fc 1 )P(fe 2 ) + 6/ NL a(fc 3 ) P [f 1 j P ^ fc2 . ) +2 eye. 
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with the terms going from A to E matching exactly Eq.(5.8) of BSS. As above, new terms appear: F 
is due to tidal term s 2 and it is generated by (s 2 c)W<5G)}, regardless of the presence of PNG, while G 
is due to n 2 , sourced by (n 2 5^5^). 
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